library(MASS)
library(sfsmisc)
snp<-read.table("SNP.txt",sep="\t",head=T)
rownames(snp)<-snp[,1]
gene<-read.table("GE.txt",sep="\t",head=T)
data2<-read.table("rlmlist.txt",sep="\t",head=T)
colnames(data2)<-c("rsid","geneid")
rownames(gene)<-gene$geneid
k1<-c()
k1<-c("gene","snp","Intercept_Estimate","snp_Estimate","Intercept_Std.Error","snp_Std.Error","Intercept_t value","snp_t value","ROBUST_F_TEST_RLMFIT")
z<-1
while(z<nrow(data2))
{
	print(z)
	kk<-c()
	flm<-rlm(as.numeric(as.character(snp[as.character(data2[z,1]),2:316])) ~ as.numeric(as.character(gene[as.character(data2[z,2]),2:390])), family=gaussian(link="log"))
	kk<-c(as.character(data2[z,2]),as.character(data2[z,1]),as.vector(coefficients(summary(flm))),f.robftest(flm)$p.value)
	k1<-rbind(k1,kk)
	z<-z+1
}
write.table(k1,"rlmlist_out.txt",sep="\t",quote=F,row.names=F,col.names=F)

#filter SITE ID 29
colnames<- colnames(snp)
k2<-c()
z<-2
while(z<length(colnames))
{
	k2<-rbind(k2,unlist(strsplit(colnames[z], "\\.")))
	z<-z+1
}
col_29<-which(k2[,2]==29)
k1<-c()
k1<-c("gene","snp","Intercept_Estimate","snp_Estimate","Intercept_Std.Error","snp_Std.Error","Intercept_t value","snp_t value","ROBUST_F_TEST_RLMFIT")
z<-1
while(z<nrow(data2))
{
	print(z)
	kk<-c()
	flm<-rlm(as.numeric(as.character(snp[as.character(data2[z,1]),-col_29])) ~ as.numeric(as.character(gene[as.character(data2[z,2]),-col_29])), family=gaussian(link="log"))
	kk<-c(as.character(data2[z,2]),as.character(data2[z,1]),as.vector(coefficients(summary(flm))),f.robftest(flm)$p.value)
	k1<-rbind(k1,kk)
	z<-z+1
}
write.table(k1,"rlmlist_out_filetered_SITE29.txt",sep="\t",quote=F,row.names=F,col.names=F)
